How to forecast product sales in the R Tidyverse with Modeltime and Prophet

Dataset:

We have 3 years of daily sales of Post-its. Those lovable note snippets that hang around your desk.

We will do the following:

1. Determine which factors influence sales.

2. Build a model to forecast daily sales.

# Install pacman if needed
if (!require("pacman")) install.packages("pacman")

# load packages
pacman::p_load(pacman,
  tidyverse, openxlsx, forecast, modeltime, parsnip, rsample, timetk, broom)

#Import data
postit <- read.xlsx("/Users/anitaowens/Documents/Marketing Analytics-Data Driven Techniques with Excel/Excel Files/Chapter 31 Excel Files/POSTITDATA.xlsx", skipEmptyRows = TRUE)

#Check results
str(postit)

'data.frame':   1096 obs. of  6 variables:
 $ Month      : num  1 1 1 1 1 1 1 1 1 1 ...
 $ Day#       : num  1 2 3 4 5 6 7 8 9 10 ...
 $ Day        : num  40544 40545 40546 40547 40548 ...
 $ Price      : num  7.52 7.52 5.95 6.2 6.1 6.2 6.98 5.95 7.12 6.98 ...
 $ Display?   : num  1 0 0 0 0 0 1 0 1 0 ...
 $ actualsales: num  390 344 636 483 486 490 524 620 416 464 ...

We have 6 variables:

1. Month
2. Day number
3. Day (This is our date variable but since we imported from Excel we get the Excel formatted dates)
4. Pricing for that particular day (7 different prices)
5. Display (This is a binary variable. 1 indicates that our product was on a display for that day. 0 indicates that it was not on display).
6. Actual sales

#Create New Formatted Date Column
postit$new_date <- as.Date(postit$Day, origin = "1899-12-30")

#Check results
str(postit$new_date)

 Date[1:1096], format: "2011-01-01" "2011-01-02" "2011-01-03" "2011-01-04" ...

#Column renaming
postit <- postit %>% 
  rename(
    day_num = 'Day#',
    display = 'Display?'
  )

#Check results
names(postit)

[1] "Month"       "day_num"     "Day"         "Price"      
[5] "display"     "actualsales" "new_date"   

#Create list of variables that need transformation
fac_vars <- c("Month", "display", "Price")

#Factor Month, Display and Price Variables
postit[,fac_vars] <- lapply(postit[,fac_vars], factor) 

#Check results
str(postit)

'data.frame':   1096 obs. of  7 variables:
 $ Month      : Factor w/ 12 levels "1","2","3","4",..: 1 1 1 1 1 1 1 1 1 1 ...
 $ day_num    : num  1 2 3 4 5 6 7 8 9 10 ...
 $ Day        : num  40544 40545 40546 40547 40548 ...
 $ Price      : Factor w/ 7 levels "5.95","6.1","6.2",..: 7 7 1 3 2 3 4 1 5 4 ...
 $ display    : Factor w/ 2 levels "0","1": 2 1 1 1 1 1 2 1 2 1 ...
 $ actualsales: num  390 344 636 483 486 490 524 620 416 464 ...
 $ new_date   : Date, format: "2011-01-01" ...

#Drop Excel formatted Day variable
postit <- postit %>% 
            select(-Day)

#From timetk package - visualize sales
postit %>%
  plot_time_series(new_date, actualsales, .interactive = TRUE)

Our plot of sales indicates that there is a positive trend of sales for our post its.

Trend in time series data is a specific pattern which in which observations are either increasing or decreasing over time. Trends can constantly change.

Task 1: What factors influence sales

#Let's do some exploratory data analysis (EDA)


ggplot(data = postit, aes(x=actualsales)) + geom_histogram()

`stat_bin()` using `bins = 30`. Pick better value with
`binwidth`.

#Sales are a bit skewed to the right

ggplot(data=postit, aes(x=log(actualsales))) + geom_histogram()

`stat_bin()` using `bins = 30`. Pick better value with
`binwidth`.

#Sales now look more normal after log transformation

ggplot(data = postit, aes(x=Month, y = actualsales)) + geom_boxplot()

#Sales are lowest in March. Sales are generally higher in Q4 and the highest in December. December also appears to have the largest spread. Some outliers during months August and September.

ggplot(data = postit, aes(x=Price, y = actualsales)) + geom_boxplot()

#This is quite interesting, as sales are highest when the price is 5.95. Sales decrease when price increases to 6.1 and so and so forth where the lowest sales happen when the price is at its highest 7.52. There appears to be a relationship between price and sales.

ggplot(data = postit, aes(x=display, y = actualsales)) + geom_boxplot()

#More sales when product is on display. 

postit %>%  
  group_by(display) %>% 
  summarize(mean_sales = mean(actualsales), median_sales = median(actualsales), sd_sales = sd(actualsales))

#Average sales are 639 vs 587 when no display.

Now we have some idea of the factors that might influence sales, let’s check for significance by using a linear regression model.

# Model Spec
model_spec_lm <- linear_reg() %>%
    set_engine("lm")

# Fit Linear Regression Model
model_fit_lm <- model_spec_lm %>%
    fit(actualsales ~ Month + Price + display, data = postit)

#Print summary of model in a tidy object
lm_summary <- tidy(model_fit_lm) %>% 
              mutate(significant = p.value <= 0.05)
lm_summary

#Use the glance function from the broom package to get additional information from the model (e.g. model statitics like r.squared)
glance(model_fit_lm)

A generic interpretation of a linear regression model. For every one unit change in coefficient (term column), moves unit sales by the value of our estimate while all other variables remain at the same level.

Most months are significant but December is positive and the most statistically significant with the highest average sales. March is the month with the lowest average sales. Time of year does impact sales.

All price coefficients are negative and significant which means that all price levels above 5.95 reduces sales and continues to reduce at each price increase. Customers seem to be price sensitive when it comes to post-its. As suggested by our EDA, there is a negative relationship between price and sales.

For display advertising, we can expect an increase of 60 when there is a display.

As for the model statistics, we have a very low p-value and a pretty decent r.square of .885, which tells us that that 89% of the variation in sales is explained by our independent/explanatory variables.

We built this linear model for mostly gaining insights into what factors influence sales.

Task 2: Build a forecast model for sales

We will split our dataset into test and training data. We will use the last 3 months of the dataset as the training set.

#Split data into test and training set
splits <- postit %>%
  time_series_split(date_var = new_date, assess = "12 months", cumulative = TRUE)

#Visualize test train split
splits %>%
  tk_time_series_cv_plan() %>%
  plot_time_series_cv_plan(new_date, actualsales, .interactive = TRUE)

We want to define our model engine. We will be using prophet.

# Model Spec
model_spec_prophet <- prophet_reg() %>%
    set_engine("prophet")

# Step 1: Model fit
model_fit_prophet <- model_spec_prophet %>%
    fit(actualsales ~ ., data = training(splits))

Disabling daily seasonality. Run prophet with daily.seasonality=TRUE to override this.

model_fit_prophet

parsnip model object

Fit time:  835ms 
PROPHET w/ Regressors Model
- growth: 'linear'
- n.changepoints: 25
- changepoint.range: 0.8
- yearly.seasonality: 'auto'
- weekly.seasonality: 'auto'
- daily.seasonality: 'auto'
- seasonality.mode: 'additive'
- changepoint.prior.scale: 0.05
- seasonality.prior.scale: 10
- holidays.prior.scale: 10
- logistic_cap: NULL
- logistic_floor: NULL
- extra_regressors: 19

#Step 2: Put model into a modeltime table
models_tbl <- modeltime_table(model_fit_prophet)

models_tbl

# Modeltime Table

#Step 3: Calibrate model
calibration_tbl <- models_tbl %>% 
  modeltime_calibrate(new_data = testing(splits))

calibration_tbl

# Modeltime Table

modeltime_accuracy() function gives all accuracy metrics.

I like to pay attention to the the mean absolute percentage error (MAPE) which is 5.95 so our forecasts are off around 6%.

#Step 4: Get Accuracy Metrics
calibration_tbl %>%
    modeltime_accuracy()

#Step 5: Create future forecast with confidence intervals
(forecast_tbl <- calibration_tbl %>%
    modeltime_forecast(
        new_data    = testing(splits),
        actual_data = postit,
        keep_data = TRUE #Include the new data (and actual data) as extra columns with the results of the model forecasts
    ))

#Step 6: Plot modeltime forecast
plot_modeltime_forecast(forecast_tbl)

Forecast sales into the future

#Create a tibble of observations with length out being the number of observations we want in reference to our date variable (new_data). Since new_date ends on Dec 31st, our future_frame starts on Jan 1st counts 90 days into the future.

#Create tibble of dates using future_frame() from timetk package
dates <- postit %>% 
  future_frame(new_date, .length_out = "1 year")

#Simulate display data
display <- rep(0:1, each = 2, length.out = 365)

#Simulate Price data
Price <- rep(c(5.95, 6.1, 6.2, 6.98, 7.12, 7.32, 7.52), length.out = 365)

#Put data into dataframe
explanatory_data <- cbind(dates, display, Price)

#Add additional Month and day_num variables
explanatory_data <- explanatory_data %>% 
  mutate(Month = format(new_date, "%m"),
         day_num = format(new_date, "%d"))

#Need to factor variables 
fac_vars <- c("Month", "display", "Price")

#Factor Month, Display and Price Variables
explanatory_data[,fac_vars] <- lapply(explanatory_data[,fac_vars], factor) 

#Change day_num from character to numeric vector and strip leading zeros
library(stringr)
explanatory_data$day_num <- as.numeric(str_remove(explanatory_data$day_num, "^0+")) 

#Do the same for the Month variable
explanatory_data$Month <- as.factor( str_remove(explanatory_data$Month, "^0+"))

# explanatory_data$day_num <- as.numeric(explanatory_data$day_num)

#Check results
str(explanatory_data)

'data.frame':   365 obs. of  5 variables:
 $ new_date: Date, format: "2014-01-01" ...
 $ display : Factor w/ 2 levels "0","1": 1 1 2 2 1 1 2 2 1 1 ...
 $ Price   : Factor w/ 7 levels "5.95","6.1","6.2",..: 1 2 3 4 5 6 7 1 2 3 ...
 $ Month   : Factor w/ 12 levels "1","10","11",..: 1 1 1 1 1 1 1 1 1 1 ...
 $ day_num : num  1 2 3 4 5 6 7 8 9 10 ...

CAUTION: It is super important that the data in your new data frame matches the exact same formatting of the data in the data used in building the forecast. If errors occur during either the model fit or forecasting phase, differently formatted data may be the culprit.

#Step 5 (Alternative): Specify the H or horizon to get a forecast into the future, but if using xregs (independent regressors), you must create a new dataframe with the xregs.

#forecast on the new tibble
forecast_tbl_future_data <- calibration_tbl %>%
    modeltime_forecast(
        new_data    = explanatory_data
    )

#plot forecast
plot_modeltime_forecast(forecast_tbl_future_data, .interactive = TRUE)

# Install ggpubr if needed and load library
if (!require("ggpubr")) install.packages("ggpubr")
library(ggpubr)

#Subset the date and prediction columns
final_fc <- forecast_tbl_future_data %>%
  select(.index, .value) %>% 
  mutate(month_year = format(.index, "%m-%Y"))

#Group by month and sum sales
final_fc_gr <- final_fc %>% 
  group_by(month_year) %>% 
  summarize(sales_forecast = sum(.value)) 

#Put forecasts into a nice and pretty table
fc_ggtable <- final_fc_gr %>% 
  ggtexttable(theme = ttheme("mBlackWhite"), rows = NULL)

#Visualize table
fc_ggtable

FINAL THOUGHTS:

I really enjoyed performing time series analysis using Modeltime package. As a Tidyverse user, this way of working with machine learning models is streamlined, functional and flexible.